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Abstract 

We present a class of non-standard numerical schemes which are modifications 
of the discrete gradient method. They preserve the energy integral exactly (up to 
the round-off error). The considered class contains locally exact discrete gradi- 
ent schemes and integrators of arbitrary high order. In numerical experiments we 
compare our integrators with some other numerical schemes, including the stan- 
dard discrete gradient method, the leap-frog scheme and a symplectic scheme of 
4th order. We study the error accumulation for very long time and the conservation 
of the energy integral. 
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1 Introduction 

Geometric numerical integration consists in preserving geometric, structural and phy- 
sical properties of the considered differential equations. Our aim is to improve the 
accuracy of geometric integrators, modifying them in an appropriate way, without los- 
ing their excellent qualitative properties (including the long-time behaviour, stability 
and the energy conservation). In this paper we focus on the discrete gradient scheme 
for one-dimensional Hamiltonian systems. 

Discrete gradient numerical schemes have been introduced many years ago in order 
to integrate numerically A^-body systems of classical mechanics with possible applica- 
tions in molecular dynamics and celestial mechanics [18] (see also lfl4l [131 [TTl l27l ). 
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Discrete gradient schemes preserve exactly (up to round-off errors) both the total en- 
ergy and angular momentum. More recently discrete gradient methods have been ex- 
tended and developed in the context of geometric numerical integration GUI . Quispel 
and his coworkers constructed numerical integrators preserving all integrals of motion 
of a given system of ordinary differential equations l2Tll22l|25ll26ll . 

In general, geometric numerical integrators are very good in preserving qualitative 
features of simulated differential equations but it is not easy to enhance their accuracy. 
Symplectic algorithms can be improved using appropriate splitting methods J2l [19l 
|28ll . Our research is concentrated on improving the efficiency of the discrete gradient 
method (which is not symplectic). 

In this paper, we continue our earlier research j9l[l0l[TT), extending the theoretical 
framework on arbitrary one-dimensional Hamiltonian systems. Numerical experiments 
are carried out in the case of the simple pendulum equation for extremaly long times 
(we test even 100 millions of periods). We study the accuracy of our new methods, 
namely, the accumulation of the global error and conservation of the energy integral. 



2 Non-standard discrete gradient schemes 

In this paper we confine ourselves to one-dimensional Hamiltonian systems 

x = H p , p = -H x , (1) 

where H — H(x,p) is a given function, subscripts denote partial differentiation and the 
dot denotes the total derivative respect to t . The Hamiltonian H(x,p) is an integral of 
motion (the energy integral). 

We consider the following class of non-standard (compare [23 1) discrete gradients 
schemes. 

x„ + i -x„ _ H (x n+ 1 , p n+ 1 ) + H (x n , p n+ 1 ) - H (x n+ 1 , p n ) - H (x„ ,p„) 
<S„ 2{p n+ \-p n ) 

(2) 

Pn+\ -p n _ H(x„,p n+ i) + H(x„,p n ) -H(x n+ i,p n+ i)-H(x n+ i,p n ) 

where 8„ is an arbitrary positive function of h,x n ,p n ,x n +i,p n +i etc. (the time step is 
denoted by h). The subscript n indicates that 8„ may depend on the step n. In the 
separable case, i.e., H = T(p) + V(x), the scheme (O becomes 

X,i+1 -x n _ T(p„+i)-T(p„) 

<5„ Pn+l—Pn 

(3) 

Pn+l ~Pn _ V(Xn+l)-V(x n ) 
8n x n +i x n 

In numerical experiments we mostly test the case T(p) = jp 2 , where further simplifi- 
cation occurs, see ( [Tot . 
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The system (O is a consistent approximation of (Q~|i if we add the condition 

5 n 

lim — = 1 . (4) 

h->0 h 

The case 8„ = h yields the standard discrete gradient method (GR), fl4l [171 [T8l l27l . 

Theorem 2.1 The numerical scheme (0 preserves the energy integral exactly ( up to 
round-off error), i.e., H(x n+ \,p n+ i) = H(x n ,p n ). 

Proof: The system l[2j implies the equality of both numerators on the right-hand sides of equa- 
tions l[2}. This, in turn, yields the theorem immediately. □ 



Therefore, any 8„ satisfying non-restrictive condition yields an energy-pre- 
serving numerical scheme. The main idea of this paper consists in finding 8 n such 
that the resulting numerical scheme is better than the standard gradient method. We 
consider and test two possibilities. First, the so called locally exact discretizations 
(section|4|, then we show that the class contains integrators of arbitrary high order. 
The corresponding /z-series for 8„ is defined in a recurrent way (section[5]l. 



3 Exact discretization 

We consider an ordinary differential equation (ODE) with a general solution x(f ) (sat- 
isfying the initial condition x(?o) = Xo), and a difference equation with the general so- 
lution x„. The difference equation is the exact discretization of the ODE if x„ = x(f„). 

It is well known that any linear ODE with constant coefficients admits the exact 
discretization in an explicit form [24], see also [1, 8, 23 1. We summarize these results 
as follows. 

Theorem 3.1 Any linear equation with constant coefficients, represented in the matrix 
form by 

^=Ax + b, (5) 
dt 

(where x = x(f) 6 W, b = const £ R" and A is a constant n xn matrix) admits the 
exact discretization given by 

«„+!=«**«„+(«**-/) A^b, (6) 
where h — t n +\ — t„ is the time step and I is the identity matrix. 

Proof. The general solution of lO is given by 

x{t) = e ('-'°' A (x(* ) +A- 1 b\ A % . 
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Taking into account that that x„ = x(t„) and, in particular, xq = x(fo), we get 

Xn = e ('n-t0)A Lu^+A^b) -A-'b , 



x„ +1 =^'»-'» + ' , ^(x(fo)+A- 1 b)-A- 1 b = e ' !/1 (x„+A- 1 b)-A- 1 b, 
which ends the proof. 



Example 3.2 Exponential growth equation: x = ax. Exact discretization: x n+ \ = e ah x n 
(a geometric series). Equivalent form: 

x n+l~ x n sn\ e al ' — 1 

= ax„ , 0(h) = . (7) 



8(h) 



8(h) 

Note that lim ——^ = 1. 

h^>o h 



Example 3.3 Harmonic oscillator: x + (Ox = 0, p = x. Exact discretization: 

x n+ i -2cos(0)/!)x„+x„_ 1 = 0, p„ = — — — ■ . (8) 

sin((Brt) 

Equivalent fo rm: 

527- + (0 z x„ + 0, 8(h ) = — sin— . (9) 

O z (rt) (B 2 

A^ofe that 8(h) w /z/or ft w 0. 

Exact discretization seems to be of limited value because, in order to apply it, we 
need to know the explicit solution of the considered system. However, there exist 
non-trivial applications of exact discretizations. In the case of the classical Kepler 
problem we succeeded to use the exact discretization of the harmonic oscillator in two 
different ways, obtaining numerical integrators preserving all integrals of motion and 
trajectories lH|5]|6l. Another fruitful direction is associated with the so called locally 
exact discretizations f5ll7l [l0ll . see the next section. 



4 Locally exact discrete gradient schemes 

First, we recall our earlier results concerning the case H = \p 2 +V(x), see JUHO). We 
tested the following class of numerical integrators 

X = t (Pn+l+Pn) ■ 

O n 2 

(10) 

Pn+\ ~ Pn _ V(x„+l)~V(x n ) 
8 n x n+l %n 
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where <5„ is a function defined by 




(id 



and, in general, x may depend on n. For simplicity, we formally assume V"(x) > 0. 
However, in the case of non-positive V"(x) one can use the same formula (either, for 
V"(x) < 0, the imaginary unit cancels, or, for V"(x) = 0, we compute the limit (0„ — > 
obtaining 8„ = h), for details and final results see iflOll . 

The simplest choice is x = xq, where V'(xo) — (small oscillations around the 
stable equilibrium). In this case 8„ does not depend on n. The resulting scheme was 
first presented in [9|, here we propose to name it MOD-GR. In [ 1 1 we considered 
the case x = x n (which will be called GR-LEX) and its symmetric (time-reversible) 
modification x — \{x n +x„ + \) (GR-SLEX). In both cases x is changed at every step. 

Definition 4.1 A numerical scheme x n+ i = l P(x„,/!)/or an autonomous equation x = 
F(x) is locally exact if its linearization around any fixed x is identical with the exact 
discretization of the differential equation linearized around x. 

We use local exactness as a criterion to select numerical schmemes of high accu- 
racy from a family of non-standard integrators, e.g. from (0. Our working algorithm 
to derive such "locally exact modifications" of numerical integrators of the form © 
assumes that 8„ depends only on x, p (or, in more general case, on x) and h. The fol- 
lowing theorem extends results of ll5l [T0ll on the case of the general time-independent 
Hamiltonian H = H(x,p). 

Theorem 4.2 The discrete gradient scheme (f2| with 



(where (O n is evaluated at x,p) is locally exact. 

Proof: We have to linearize the continuous system Q, then to find the exact discretization of the 
obtained linearization. Therefore, we put x = x + %, p = p+r) into Q} and neglect all terms of 
order greater than 2. Thus we get 





(12) 



£=H P +H px $+H pp ri , 



f] — —H x —H xx t; —HxpTj . 



(13) 



The exact discretization of the system d!3t is given by 




(14) 



(compare Theorem l3.lt . where 




(15) 



5 



We proceed to the linearization of the discrete system {2). We substitute 

X„ = X + £,„ , p n =p+ T] n , (16) 

and assume that 8„ depends only on x, p and h (which is equivalent to taking only the first, 
constant, term of the Taylor expansion of 8„ with respect to ^ n ,T)„). Then, we linearize the 
system ([2} around x, p (neglecting terms of at least the second order with respect to £,„ and rj„), 
obtaining 

— — = H p + -H xp (|„ + 1„ +] ) + -jH P p (r/„ + r/„ + i ) , 

1 1 

" + g = —H x - -H xx (£„ + |„ + i ) - -H xp (r/„ + r/„ +1 ) , 

where partial derivatives H x ,H p ,H xx ,H xp and H pp are evaluated at x,p. After simple algebraic 
manipulations we rewrite this linear system in the form 

6h-i ^ =m ( 6. )+w _ (l8) 



where 



1 + I co 2 S 2 ^ - S„//„ 1 - S„H xp - \ to 2 S 2 



8„ ( \ + \8 n H xp \8 n H pp ^ I H p 



(19) 



l+l^n^n \ —\$nH xx 1 — j8 n H xp J \ —H x 

and co„ is defined by d!2t . Taking into account J 1 5b . we get 

l + ±o 2 S 2 1+K^ 2 ' 1 + W 5 « 2 1 + Ja#j? 

Systems dl4b and J 1 8b coincide if and only if 

M = e bA , v=U A -i}A~ l h. (21) 

The proof reduces to showing that the system d21t is identically satisfied if 8„ is given by jilt . 
We easily verify that 

A 2 = -ffl 2 /, al = H xx H pp -Hl p . (22) 

Hence 

e hA = cos hco„ + a>~ l A sin hco„ . (23) 

The second equation of d2 1 b is satisfied for 8„ of any form. Indeed, using the first equation of 
d21t and then the first equation of {22}, we get 

(M-/A~ 1 b= V 2 , " " ' = } " ' =w. (24) 
1 + > 2 S 2 1 + 

Finally, the first equation of d2 1 b is satisfied if and only if 

!-W 5 k u S " sinhco n 

-cosAffl„, = — — - = , (25) 



1 + ' l + ico 2 5 2 co„ 
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where we took into account d20t and d23 b . From the first equation we compute 

1 1+ COS ft©,, 2 ll( °n 

— 1 ,-, = 5 = cos — , (26) 

and substituting it into the second equation of d25t we get J 1 2b - □ 

Remark 4.3 Assuming x = x n , p = p„ we get a numerical scheme called GR-LEX, 
while the choice x = j(x n + x n+ \), p = j (p n + Pn+i) yields another scheme, named 
GR-SLEX. The system (0) is symmetric (time-reversible). The numerical scheme GR- 
SLEX preserves this property, while GR-LEX does not preserve it. 

The discrete gradient schemes GR and MOD-GR are of second order. Locally exact 
discrete gradient schemes have higher order: GR-LEX is of 3rd order and GR-SLEX 
is of 4th order, see flOl . In the next section we show how to construct discrete gradient 
schemes of any order. 



5 Discrete gradient schemes of Mh order 

We consider the family (|2]) of non-standard discrete gradient schemes for the Hamilto- 
nian system ((T). The family is parameterized by a single function 8„ and this function 
can be expressed by x„,p n ,x n+ \,p„ + { as follows, 

8 = 2Q„ +1 -x n ){p n+ i -p n ) 

H(x„ + i,p n+ \)+H(Xn,Pn+l)-H{x n+ \,p„)-H(Xn,Pn) 

Replacing here x n+ \ : p n+ \ by the exact solution x(t n+ i),p(t n+ \) we formally obtain 
8 n corresponding to the exact integrator. In practice, we can replace x„ + \,p n+ i by 
truncated Taylor expansions (and truncate the final result). 

Therefore, we take Taylor expansions truncated by neglecting terms of order higher 
than jV (see Appendix A, formulae d3T1 >) and compute 

2(x [N] -x )(d [N] -d) n 

ff(4+i,pKi)+^(^/'Ki)-^Si ) /'«)-^(^. J pn) ^ 

where coefficients are functions of x n and p n . Then, truncating the obtained result, 
we define 

8f ] = f j a k (x n , Pn )h k . (28) 

k=\ 

The first few coefficients reads 

a\ = 1 , «2 = , «3 = H XX H PP — H%p — H x H xp p — H p H xxp , 
a 4 = H^H xpP p — H^Hxxxp + HpH p pH xxx — H x H xx Hp PP — 3H p H xp H XX p (29) 
~~^3H x H X pH X pp , 
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Figure 1: Energy as a function of time (t = nh), h = 0.25, po = 1.8, E ex = 0.62. 

where all partial derivatives are evaluated at x = x„, p = p n . In the separable case, 
H = T(p) + V(x), the formulae simplify 

a\ = 1 , ai = , «3 = j^TppVja , a/\ — j^ {T p T pp V^, x - V x V xx Ts p ) , 

a 5 = (9V 2 V xx T Ap + 9T*T pp V Ax - 12V x V 3x T 2 p - 12T p T 3p V 2 (30) 

+ 6^2-16^,^3.0 

where Vk x denotes kth derivative of V with respect to x, etc. The case H — jp 2 + V (x) 

is discussed in more detail in ifTTl . where explicit formulae for Sn f° r N ^ 1 1 can be 
found. 

The gradient scheme (O with 8„ = is called GR-N. Its order is at least N, 
sometimes higher (e.g., GR-1 is of 2nd order). Actually GR-1 and GR-2 are identical 
with GR. 



6 Numerical experiments 

The accuracy of high-order discrete gradient schemes was tested on the case of the 
simple pendulum, H = jp 2 — cosx (for simplicity always assuming xq = 0). We com- 
pared GR-3, GR-7 and GR-LEX with the discrete gradient method (GR), the leap-frog 
scheme (LF), 4th order explicit Runge-Kutta method (RK-4), high-order Taylor meth- 
ods (TAY-N, see appendix A) and a 4th order symplectic scheme (SP-4, see appendix 
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Figure 2: Energy as a function of time (f = nh), h = 0.25, po = 1.8, E ex — 0.62. The 
line E = 0.62 corresponds to GR, GR-LEX, GR-3 and GR-7. Other, scattered, points 
are produced by TAY-10. 

B). Computing global errors we use the exact solution of the simple pendulum equa- 
tion, expressed in terms of elliptic integrals. 

In previous papers |9 10, 1 1 1 we focused on the stability and accuracy of the period 
(all motions of the pendulum are periodic). Here, we test the global error, accumulated 
after 120 periods (Figures|4]and|5]) and the accumulation of error after a very long time 
(up to n = 10 8 steps), Figure [6] We also check the preservation of the energy integral 
by different numerical schemes, Figures Q] [2] and [3] Details concerning the solution 
of implicit equations are the same as in [10], e.g., at every step we iterated until the 
accuracy 10~ 16 was obtained. 

Symplectic methods are known to preserve almost exactly the energy integral |fl5l , 
some positive results in non-symplectic case are also known [ 12|. FigureQ]shows how 
accurate is the preservation of the energy by symplectic integrator SP-4 as compared 
with TAY-5 (permanent, fast growth of the energy) and RK-4 (energy is decreasing ap- 
proaching the stable equilibrium value). A high order of a given scheme is not sufficient 
to assure the conservation the energy. From the beginning TAY-10 produces small, but 
permanent, drift of the energy, while all discrete gradient schemes yield almost exact 
value of the energy, see Figure[2] According to Theorem l2.1l all gradient schemes pre- 
serve the energy integral exactly (up to round-off errors). Only after very long time 
one can notice that also discrete gradient schemes have a slight drift of the energy. A 
curious phenomenon can be observed at Figure[3] Gradient schemes and TAY-10 show 
a linear growth of the energy error (but the energy error of TAY-10 is always greater 
by 6 orders of magnitude!), while the energy errors of symplectic schemes, LF and 
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Figure 5: Global error as a function of h, evaluated at t = 1207}/, for po = 1.8. 




SP-4, vary in a large range but do not show any systematic time dependence. However, 
in the considered time interval, discrete gradient methods preserve the energy more 
accurately by several orders of magnitude than symplectic integrators like LF or SP-4. 

Figures[4][5] show /t-dependence of the global error calculated at t = \2QT ex , where 
T ex is the period of the exact solution (e.g., T ex — 6,283342396 for po = 0.02 and 
Tex = 9.12219655 for po = 1.8). We observe that, usually, higher order integrators are 
more accurate. An important exception is GR-LEX, of 3rd order, which for po = 0.02 
and h > 0.3 is better than GR-7 (the behaviour of GR-SLEX is almost the same as 
GR-LEX). However, TAY-10 is clearly the best in this case. Only after very long time 
evolution TAY-10 becomes less accurate than GR and SP-4, although initially it was 
comparable with gradient schemes of high order, see Figure [6] Several schemes at 
Figure [6] show linear growth (at least for large t). It has been shown, see [3] [15], that 
symplectic integrators (under some mild conditions) have linear error growth. Results 
of our experiments suggest that after sufficiently large time some other schemes (e.g., 
RK-4, GR, GR-3) also accumulate error linearly. Finally, we point out that until n = 10 7 
the global error of GR-LEX is smaller than the period (and the error of GR-7 is even 
smaller, by one order of magintude). The global error of GR is smaller than that of 
SP-4, not saying about GR-3 or GR-7, see Figure[6] 

7 Conclusions 

Modifications presented in this paper essentially improve the discrete gradient method 
(in the one-dimensional case) keeping all its advantages. Modified gradient schemes 
GR-LEX, GR-SLEX, GR-N have important advantages: 

• conservation of the energy integral (up to round-off errors), 

• high stability, exact trajectories in the phase space, 

• high accuracy (third, fourth and Mh order, respectively), 

• very good long-time behaviour of numerical solutions. 

We point out, however, that numerical schemes like all discrete gradient methods, 
are neither symplectic nor volume-preserving. Most of them, including GR-LEX and 
GR-N (N > 2) are not symmetric (time-reversible). GR and GR-SLEX are symmetric. 

In the near future we plan to generalize the approach presented in this paper on 
some multidimensional cases [7| (the crucial point is that 8„ is a matrix) and to extend 
the range of its applications on some other numerical integrators (including the implicit 
midpoint rule and numerical schemes which preserve integrals of motion [21 25]). One 
can also use a variable time step, if needed (7). 
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Appendix A. Explicit Taylor schemes of Mh order 

Mh Taylor method for the system (Q~|i is defined by 



» b k h k 

X n+l = L 



k=0 



k\ ' 



[N] _ y c k h k 

Pn+1 — 2-i ft ' 



k=0 



(31) 



where the coefficients b k , c k are computed from Taylor's expansion of the exact solution 
(see, for instance, 1161 . p. 18). We assume x(t) — x n , p(t) — p„, t — t„ and expand 
x(t + h) and p(t + h) in Taylor series: 



Uj.u\ V **** 



k=Q 



k=0 



h = 



Ck 



d k x(t) 



dt k 



d k p{t) 



dt k 



(32) 



(33) 



where all derivatives are replaced by functions of x n ,p„ using ([T} and its differential 
consequences, e.g., 



d dH 

X = — ^r— = H px x + Hppp — H px H p — HppH x 



dt dp 

where H x ,H p ,H xx etc. are evaluated at x n ,p n . Thus we get 



(34) 



h=H p , 



bi — H p H xp - H x H pp . 



b3 — H^H 3p + H xxp Hp — 2H x HpH ppx + H p H xp - H p H pp H 



c\ = -H x 



ci — H x H xp — H p H x 



(35) 



C3 = —HpHix — H xpp H x + 2H x H p H xxp — H x H^ p + H X H PP H XX , 
and subsequent coefficients can be easily computed using the total derivative: 
db k db k db k dc k dc k 

bk+l = ~dT = p ~dx~~ x lfp ■ Ck+l =Hp ^x~- Hx Jp~ ■ 



(36) 



The formulae d35l ) simplify in the case H = T{p) + V(x) fTD . 
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Appendix B. Explicit symplectic schemes of 2Mth order 

Symplectic explicit integrators of arbitrary even order N = 2M can be derived by com- 
position methods, see, e.g., llT5l . In this section we present results of the pioneering 
paper ll28l . confining ourselves to the case H = jp 2 + V (x). 

The numerical scheme SP-2M is defined by the following procedure. Having x n , p n 
we compute the next step, x n+ \,p„ + \, as follows. We denote x„ =x™,p n — pt°] and 
perform K+ 1 iterations (where K = 3 M ~ ) 

xW = jtt 1 - 1 ! + hc^p^ , pW = p^ - hd l £V(j<M) , (37) 



where cl(£ (i = 1,2, . . . ,K+ 1) have to be carefully computed (see below) in order 
to secure the required order. Then we identify =x„ + i,p^ +1 l = p„+\. 

The coefficients cj^,^ are computed recursively. First, all coefficients for M = 1 
are given by: 

4 1] = ~, 4 1] = i, 4 21 = 5, 4 ] =o. 08) 

Then, we express coefficients c^ +1 , d® +1 by coefficients cj,'J , d m '■ 



d® ~ 

"m+1 — 


[2*4 
"m+1 


'] - v M a - 1 

— }m a m , {' — 1 j 


...,*), 


Ak+i] _ 
"m+1 ~ 






.,*), 


j[3*+l] 
"m+1 


= o, 






c M - 

L m+1 — 


J 2 *+ 
c m+l 


<+i] _ [i+i] f . . 


= 1,2,...,*:), 


L m+\ 


= (1 


-2y m )c|,' ! +11 , (i = 


1,2,...,*-!), 


_ 

L m+1 ~ 


= (1- 


v )(c ll] +ct 1] ) 




m+1 


= (1- 


-v )(c [l] +c [k+l] ) 

jm) 1 c m L m J j 





2(2-2V3) ' 2 2 2(2- 2^3) ' 



*4 ' ~ *A ' — 2 _2i/3 ' *4'~~ 2-2 1 / 3 ' *4' — 0- 



(39) 



where k = 3 m 1 and ^ 

- Vm = 2-2 1 /(2m+l) ' ^ 40 ) 

In particular, the symplectic integrator SP-4 has the following coefficients 

Ji] _ r W _ 1 M _ M _ 1 - 2 ' /3 

2 — c 2 — t/o t1/Ti > L 2 — c 2 



(41) 



The scheme SP-4 was independently presented in [ 1 3 1 and |28|. 
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